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ABSTRACT 


The process of hydrogen transformations in 
metal -hydrogen systems has been simulated using the finite element 
method (FEM) . The aim of the finite element evaluation was to 
calculate the accommodation energies as a function of progress of 
hydrogen transformations. The a-/3 transformation in the 
niobium- hydrogen system has been specifically addressed. The 
accommodation energies during the phase transformation has been 
evaluated by assuming both linear elastic and elasto-plastic 
behaviour of the system. In elasto-plastic analysis NbH is assumed 
to be elastic and Nb is assumed to be undergoing plastic 
deformation following the linear strain hardening law. In the 
finite system considered, the transformation progress from the 
centre of the system outwards and from the surface of the system 
inwards have been addressed. Accommodation energies have been 
evaluated for both a->/3 (Nb-»NbH) and /3-»a (NbH-»Nb) transformations 
using plane stress and axisymmetric formulations. The effect of 
changing the aspect ratio of a disk shaped and prolate spheroidal 
precipitate on the accommodation energy has been investigated. 
The results obtained in elastic analysis indicates that the total 
accommodation energy decreases with increasing aspect ratio, for 
all the transformation cases considered. Moreover, most of the 
energy is accommodated in the outer phase which surrounds the 
inner phase in the elastic analyses. The accommodation energies 
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obtained in the plane stress case are lower than that in the 
axi symmetric case because of the constraint against the free 
expansion is lower. The finite element elasto-plastic analysis of 
the hydrogen transformations revealed that plastic deformation of 
Nb results in lowering of the elastic energy of the system. The 
variation of this plastic energy for all the transformations has 
been discussed as a function of the aspect ratio. 
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CHAPTER 1 
INTRODUCTION 


1.1 GENERAL INTRODUCTION 

Hydrogen transformation refers to the class of 
transformations in which the new phase precipitates essentially 
due to the diffusion of only hydrogen [1] . Hydrogen transformation 
can be further classified in to two broad groups. In forward 
hydrogen transformation, the hydride precipitates from the 
metallic matrix. This is generally observed in metal -hydrogen 
systems at low temperatures where the metal atoms in the lattice 
are relatively immobile and hydrogen diffuses with ease in the 
metallic matrix. The subsequent ordering of these mobile hydrogen 
atoms on certain fixed subsets of interstitial sites in the parent 
matrix leads to the formation of hydride phase in the matrix. 
Reverse hydrogen transformation, i.e. involving hydride 
decomposition can be visualized to occur when hydrogen is removed 
from the hydride phase which leads either to the transformation of 
a hydride to a hydrogen free structure or to the transformation of 
hydride of higher hydrogen concentration to a hydride of lower' 
hydrogen concentration. 

In metal -hydrogen systems, the precipitation of a solid 
second phase from a solid matrix results in the creation of strain 
energy in the system [2] . This strain energy arises because the 
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transformation strains cannot be relaxed in the case of solid 
state phase transformations. Transformation strains are present in 
solid state transformations due to molar volume differences 
between the precipitating phase and the matrix. The 
transformation strain energy or accommodation energy controls 
several features of the transformation. The kinetics of 
transformation and precipitate morphological characteristics are 
governed by its magnitude [3] . The free energy changes of solid 
state transformations are affected by the strain energy 
contribution [3] . Repeated formation and decomposition of 
hydrides in metal -hydrogen energy storage systems result in 
hysteresis losses [4] . This hysteresis loss has to be minimized 
while developing commercial hydrogen storage materials, because it 
leads to a loss in useful available energy. Moreover, 
accommodation energies determine the embrittling effects of 
hydrogen in many commercially important nuclear materials [5] . 

There is a volume increase upon precipitation of hydride 
from a metallic matrix, which has been observed in the Nb-H 
[6] ,V-H [7]., Zr-H [8] and Ti-H [9] systems. In metal -hydrogen 
systems, the molar volume of the hydride phase is generally 20 to 
25% greater than that of the matrix from which it forms. 
Therefore, during hydride formation, the free expansion of the 
hydride is restricted by the presence of the surrounding matrix. 
This creates compressive normal stresses at the precipitate-matrix 
interface. In a reverse hydrogen transformation, tensile normal 
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stresses develop. The presence of these stresses during the 
forward and reverse hydrogen transformations creates strain energy 
in the system. For the purely elastic case, neglecting the strain 
energy due to hydrogen solution, which occurs before the beginning 
and after the end of hydrogen transformations, the system is 
stress-free at the beginning and end of the hydrogen 
transformation. The energy associated with the volume change upon 
transformation is stored as elastic energy in the matrix and the 
precipitate. If the metal or hydride phase yields during the 
transformation, the strain energies are relaxed to lower values 
due to irreversible plastic deformation [1] . 

The transformation strain energy has been previously 
related to elastic [10] , plastic [11] and elasto-plastic [1] 
strain energy changes over the entire cycle of the hydrogen 
transformation. Most of the earlier strain energy calculations 
during solid state transformations have been performed assuming an 
infinite system size under isotropic continuum conditions. An 
infinite continuum is rarely encountered in practice. Moreover, 
isotropic conditions rarely prevail in the matrix or the 
precipitate. In the case of practical metal -hydrogen systems, the 
finite particle size leads to relaxation of energies due to the 
presence of free external surfaces. No satisfactory theory exists 
to model this effect; although, the effect of the surfaces has 
been conventionally treated using an extra 'image' stress field 
[3] . Recently, the finite element method has been used to model 
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the phase transformation in metal -hydrogen system, in which the 
effect of the free surface is also considered. In a pioneering 
work, Sen et al developed an elastic [12,13] and elasto-plastic 
[14] finite element model to calculate the accommodation energies 
during forward and reverse hydrogen transformations as a function 
of progress of transformation in the Nb-H system. The precipitate 
was considered as both a solid and hollow sphere. However 
hydrides in general tend to have ellipsoidal or plate-like shapes 
rather than spherical shape [15] . 

1.2 OBJECTIVE OF THE WORK 

The main objective of the present work is to develop a 
finite element model for calculating the accommodation energies 
during solid state hydrogen transformation, where the shape of the 
precipitate is a general ellipsoid. The finite element method is 
selected for the calculation because it is a versatile method and 
it is possible to predict, with sufficient accuracy, the 
distribution of stress, strain etc. for a complicated problem. It 
is intended to prepare a displacement -based finite element program 
employing the initial strain method. It is proposed to carry out a 
comprehensive investigation of the effect of changing the aspect 
ratio on the accommodation energies of the precipitate and the 
matrix. Both elastic and elasto-plastic analysis are undertaken 
for the plane stress as well as axisymmetric cases. The hydride 
transformations Nb-NbH and Pd-PdH have been addressed in this 
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thesis . 


1.3 THESIS ORGANIZATION 

Chapter 2 of this thesis contains the literature review. 
Theoretical solid- solid transformation models proposed by earlier 
workers have been discussed here. The concepts of the finite 
flement method (FEM) , and the models developed using FEM have been 
discussed in detail in Chapter 3 . The results obtained using the 
FEM models and their comparison with other mathematical models and 
experimental observations have been presented and discussed in 
chapter 4 . The conclusions are drawn and the scope for further 
study are presented in Chapter 5. The Appendix A provides the 
results for the Pd-H system. 
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CHAPTER 2 

LITERATURE REVIEW 


2.1 SOLID STATE PHASE TRANSFORMATION 

2.1.1 Introduction 

The free energy change (AG°) of phase transformation is 
a balance between the reduction in the volume free energy of the 
transformation and the increase in surface free energy due to the 
precipitation of a new phase during the transformation. If the 
transformation takes place in a stress free condition, and the 
nucleation of a spherical precipitate is considered, the free 
energy change of the transformation is given by [16] 

AG° = — n a^ AG + 4 n a^ Ao- . .2.1 

3 V 

where, AG^ is the the volume free energy change associated with 
nucleation free energy change per unit volume. Act is the free' 
energy change per unit surface area created and a is the radius of 
the spherical precipitate nucleus. The variation of the two energy 
terms with the size of the precipitate (a) are shown schematically 
in Figure 1 [16] . The above expression holds good for phase 
transformations involving at least one non-rigid phase so that it 
takes place under unconstrained conditions. However, in phase 
transformations involving the precipitation of a solid phase in a 
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solid matrix, the transformation involves development of stresses 
within the system [3] because the specific molar volumes of the 
matrix and the precipitate (which forms from it) are normally 
different. The free energy expression then needs modification to 
account for the transformation stresses. 

In solid state transformation, the strain energy term 
(W) has to be incorporated in the relation expressing free energy. 
The transformation free energy {AG°) is then expressed as [17] 

AG° = -^na^/3 [AG^+W] +na^A(r [2 +g(/3)] ..2.2 


where , W is the elastic strain energy per unit volume, 0 is the 
aspect ratio of the precipitate, Act is the nucleus/matrix 
interfacial energy, and the function gO) is given by hyperbolic 
functions for ^<1 or /3>1. For /3=1, g(|3)=2. For the formation of 
the critical nucleus, the required conditions are 


[5(AG°)/Sa]^ = 0 ..2.3 

and 

[5(AG°)/5^]^ = 0 . .2.4 

a. 


From equations 2, 3 and 4, the critical nucleus size (a*) for 
stable nucleus formation is given as 


cr [2 + g(|3) ] 
2 /3 [AG + W] 


From equations 2.2 and 2.5, the activation energy for the 
. * 

formation of a stable embryo (AG ) for a spherical nucleus is 
given by 
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free energy cl 



radius a 


Figure 1 . Variation of free energy of transformation as a function of 
without the strain energy contribution. 


nucleus radius 



.. 2.6 


AG* = n o-^ [2 + g(g)]^ 
12 [AG + W ] ^ 

V 


or, 


* , ■* 
AG /AG 

h 


[2 + g(p)]^ 

[8/3 (1 + W/AG^] ^ 


..2.7 


^ o o 

where, AG^ is (16 H cr /3 AG ) , the free energy for the formation 
of a spherical nucleus in the absence of strain energy. The above 
equations have been formulated for homogeneous nucleation and have 
to be slightly modified to account for heterogeneous nucleation. 


2.1.2 Continuum mechanics models 

Many models have been earlier evolved to evaluate the 
stresses, strains and strain energies involved during solid state 
transformations. They are discussed briefly below. 

2. 1.2.1 Elastic models. 

Most of the earlier models have assumed a purely elastic 
continuum. The first study of the elastic strain problem were 
conducted by Eshelby [2] whose work has since become a classic. 
The elastic strain energy induced by the coherent ellipsoidal 
inclusion in an elastic isotropic medium was calculated on the 
assumption that both the phases have the same elastic moduli. The 
general theory of strain energy was proposed by Khachaturyan [18] , 
in which The exact equation for strain energy and the fourier 
transform of the elastic displacements in an arbitrary two-phase 
coherent matrix were derived. Sass, Mura and Cohen [19] expressed 
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the displacement field associated with a cuboidal precipitate in 
the form of a fourier series. Favier [20] proposed the exact 
solution for inclusions having the shape of rectangular 
parallelepiped. Sankaran and Laird [21] calculated the deformation 
field associated with a square inclusion. All these calculations 
have been carried out for elastic strain field in the limiting 
case of isotropic elasticity and homogeneous moduli. Lee and 
Johnson [22] evaluated the strain energy of a coherent cuboidal 
precipitate in an anisotropic matrix. In the following 
paragraphs, the Eshelby method as elucidated by Khachaturyan [3] 
would be presented. 

The formation of a multi phase alloy has been considered 
and the following assumptions have been made.* 

1. The parent and transformed phase particles have the 
same elastic moduli (the homogeneous moduli approximation) . 

2 . The average length dimension of the precipitate and 
the average distance between the nearest precipitates are small 
compared with the typical length dimensions of the crystal . 

3 . External boundaries of the multiphase mixture are 
stress free. 

The stress -free strain tensors 


e 


0 

1 j 


( 1 ) , 




(p) , 



.. 2.8 


describe the macroscopic shape deformation of the parent phase 
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caused by crystal lattice re-arrangements of the 2'"'^, 

kinds associated with the phase transformation in the 
stress-free state. 

The transformation to the multiphase mixture from the 
parent matrix is modelled to involve six successive 
steps (Figure 2) . 

Step 1 : A small group of clusters is chosen in the 
parent phase crystal, and cut out of it. If the clusters are large 
enough to ignore the surface effects, the energy of the system 
remains unaffected in this step. 

Step 2 : Each cluster is allowed to transform to a new 
phase corresponding to one of the crystal lattice re-arrangements 
under stress-free conditions. The clusters are thus transformed 
into a set of new phase precipitates. The associated expansion of 
clusters involves no strain energy change since clusters are 
stress-free. The expansion is described by strain tensors c°j(p) 
(p =1, 2 ,...., v) associated with the macroscopic shape 
deformation in the stress-free state. 

Step 3 : Surface traction is applied to each cluster to 
restore the shape it had before the transformation. The 
restoration of the shape of the p’^*' kind particle requires the 
homogeneous strain 



cut off 


(step 1) 



phase transformation 
in stress-free state 
(step 2) 



introoduction 
back in to holes 
and ’welding’ 


surface traction and 
shape restoration 
(step 3) 

(step 4 and 5) 



relaxation after taking 
external forces off 


Figure 2. Successive steps to form coherent new phase inclusion [3]. 
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According to Hooke's law the internal stress 
traction is 


O' 


1 J 


= X (■ 

1 j k 1 ' 


e;/p)) 


induced by surface 

.. 2.10 


where is the elastic modulus tensor. Since the final 
homogeneous strain within a kind particle is determined by 
equation 9/ the mechanical energy changes by 


AE 


self 


(P) 


£ 

2 


V X 

P 1 Jkl 


e (p) 

1 j ^ 



.. 2.11 


where is the volume of a particle- of the p*'*' kind. The total 
energy change of all kind particles is 


N 

p 


AE 


self 


(P) 


1_ 

2 


N V 
p p 


X 


1 J k 1 





(p) 


V X 
P I Jk 1 


c (p) 

1 j ^ 


,(p) 


.. 2.12 


where N is the number of all p*^^ kind particles, V =N v is the 

P 0 - 4 . P P P 

total volume of the p kind particle. Finally, the total energy 

t h 

change associated with the shape restoration of all particles is ' 
the sum of .the energy changes given by equation 2.12," and is 
written as 


AE 

3 


1 

2 


S V 

p P 


X 


1 J k 1 


= u<P> 



..2.13 


Step 4 : The particles (p= 1, 2,...,i/) are re-introduced 
back into their holes left in the parent crystal after their 
removal. Since after step 3 each particle just fits into the space 
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from which it was removed, the insertion does not give rise to any 
energy change: 

AE = 0 ..2.14 

4 

step 5 : The particles are now welded to the parent 
phase (matrix) . 

Step 6 : Finally, the particles relax by introducing the 
equilibrium elastic strain into the lattice. During relaxation, 
each particle will initiate crystal lattice displacements in its 
vicinity. The displacement is opposed by elastic resistance of the 
lattices. The associated relaxation energy per unit volume 
J (r) may therefore be presented as power series in the 

deformations e^^(r) . Truncation after the square term yields 


f , (r) 

relax 


-O' (r) 

1 J 


1 J 


1 J k 1 


1 J 


k 1 


..2.15 


where cr° (r) and X, are the first- and second-order expansion 
1 J LlkJ ^ 

coefficients. The total relaxation energy is given as integral of 
equation 2 . 15 over the system volume V 


AE 


relax ’ 


[-cr (r)e + — X G e ]dV ..2.16 

IJ 2 IJkl IJ kl 


Since the elastic relaxation is a spontaneous process, 


AE so 

relax 


..2.17 


Physically, the effect of relaxation is to remove a part of the 
elastic distortion energy introduced in step 3 when transformed 
clusters were deformed to fit properly into the parent phase body. 
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Summation of the contributions to strain energy from all 
steps gives the total elastic strain energy: 


E 


elastic 


AE + AE 

3 relax 


I 


Z V X 

p p 


(p) (p) + 

i j k 1 i j k i 




i J 


I J k 1 


i J 




. . 2.18 


Linear terms in appear in equations 2.15, 2.16 and 

2.18 because a multi -connected heterogeneous medium is being 

considered. Unlike the usual case of homogeneous media, the 

systems under consideration are strained in the stress-free state. 

In other words, if cr (r) s o, we have c (r) * 0 (cr is the 

1 J i J i J 

stress tensor) . 

Let the strain tensor of a stress- free state 
(o"^ J (r) = 0) be 8°^ (r) . Using the notation of equation'^ we have 

e°^(l) if r is inside any first-type particle 
c^j(r) = e°j(p) if r is inside any p’^*' type particle 

if r is inside any type particle 

■0 otherwise ..2.19 


The condition 2.19 may be written in the condensed form: 

0 ° (r) = Z e (r) 0° , (p) . .2.20 

'J p = i 'J 

where 0 (r) is the shape function of particles of type p equal to 
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unity inside a particle and zero outside it. The shape function 
(r) is in general a multi connected function and may describe an 
arbitrary set of particles of the p’^^ type. 

According to the usual relation of elasticity, the 
elastic stress is related to the elastic strain as 


cr,j(r) = 


S 'P elastic 

S'Cj j (r) 


-<T ir) + X c ..2.21 

1 J 1 J k 1 k 1 


The strain, e°j(r), corresponding to the stress-free state 
may be found from equation 2.21 by setting crjj(r) a 0 


-cr (r) + X c = 0 

1 J 1 J k 1 k 1 


. . 2.22 


Substituting equation 2.20 into equation 2.22 yields the 


definition of the function cr° (r) : 

cr; (r) = S, 0 -; (p) e (r) 

1 j p = 1 1 J P 

where 


..2.23 


cr° (p) = \ c° (p) ..2.24 

1 J 1 J k 1 1 J 

Substituting equation 23')into equation 1^ gives 


AE 


relax 


V 

[ Z cr° (p) e (r) c + — A c e ] dV ..2.25 

p=l p ' ij 2 IjkllJ kl-' 


Solving these expressions, the elastic energy can be 
evaluated. 

On the lines of the Eshelby model, expression for the 
elastic strain energy of coherent ellipsoidal precipitates in 
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anisotropic crystalline solids have been derived by Lee et al. 
[8] . The model has been used to evaluate strain energies of 
isotropic precipitates in an anisotropic matrix, and anisotropic 
precipitates in an isotropic matrix. One of the assumptions of the 
theory of Eshelby have been used, namely, when an ellipsoidal 
precipitate with a constant stress -free transformation strain 
is embedded in an isotropic matrix, the constrained strain, is 

also constant inside the precipitate. The suitability of the model 
to the anisotropic case has been illustrated in the literature, 
through examples. The expression for the elastic strain energy per 
unit volume of the precipitate in a homogeneous system, has been 
derived as 


W 






1 J k 1 



.. 2.26 


In the case of an inhomogeneous system, in which the 
elastic constants of the precipitate and the matrix are not the 
same, an 'equivalent' stress-free transformation strain c°^ is 
first obtained by solving 

.. 2.27 


\ (e 

i J k 1 k 1 


e ) 

k 1 


X {c 

1 J k 1 k 1 


C V 

C ) 
k 1 


where X* is the elastic constant of the precipitate and ^ is 
the actual stress-free transformation strain of the precipitate. 
Hence, the strain energy per unit volume of such a precipitate is 
given by 


W 


L 

2 




* 

X 

1 J k 1 



.. 2.28 
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The model has been validated with data for various 
inhomogeneous (different properties of matrix and precipitate) , 
isotropic systems . The results have been reported in a normalized 
form W/W^ , where is the elastic strain energy per unit volume 
of precipitate in the isotropic homogeneous case, given by 


W 


2 M 


1 + 1 / 
1 - 1 / 


2 


C 


..2.29 


where 4 and 1 / are the shear modulus (=C ) and Poisson's ratio 

4 4 

{=C^2/2(^i 2 matrix phase. The model has been used 

to analyze the elastic strain energy variation for various 
combination of matrix and precipitate properties like shear 
strength ratios between the precipitate and the matrix, aspect 
ratio of the precipitate, and anisotropy of the precipitate and 
the matrix properties. The analysis indicates that in systems 
having high degree of anisotropy, markedly different results are 
obtained from those obtained from isotropic elasticity. 

2. 1.2. 2 Elasto-plastic models 

Recently, Lee et al. [23,24], have proposed models to 
represent solid-solid phase transformation, even though they still 
employed the basic Eshelby approach. The novelty of their models 
is that the concept of plastic deformation has been included in 
these models. This was lacking in all the earlier models. The 
models have been developed for three specific cases: the linear 
elastic case, for a matrix undergoing ideal plastic deformation, 
and for a strain hardening system with linear hardening and power 
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law hardening. These models are discussed below. 

In the pure elastic case [23] , a misfitting spherical 
precipitate, whose radius is a(l+c) in the absence of constraints, 
is introduced into a spherical hole of radius a in an infinite 
matrix with shear modulus u and Poisson's ratio v . The precipitate 
centre is taken at the centre of the system. The precipitate is 
assumed to change to a(l+^e ) under the constraint of the matrix. 
Assuming spherical symmetry and the equilibrium equations of 
elasticity, the strain energy per unit volume of precipitate (note 
that as per the convention of the model, the volume of the 
precipitate, that has been used to calculate the energy per unit 
volume of system, is equal to the volume of the matrix from which 
it forms) is given by 


E 


e 1 


6 4 g y 

g ( y - 1) + 1 


..2.30 


where y is the ratio of the bulk moduli of the precipitate and the 
matrix, p is the shear modulus of the matrix, e is the stress-free 
misfit strain and g = (l+i^) / (3 (1-v) ) . 

Plastically deforming systems have been modelled with 
the assumption that the precipitate is incapable of plastic 
deformation under the hydrostatic stresses acting on it [23] . The 
plastic deformation takes place in the surrounding matrix (Figure 
3) . In the model, a and r are taken to be the precipitate radius 

P 

and the plastic zone radius, respectively. The Von mises yield 
criterion has been utilized. The total elasto-plastic strain 
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energy for a system exhibiting ideal plastic behavior is given by 
[23] 


w = 


^2 O' 

P , y 


2rK 


(i-r) + 


an 


( 1 

f r N 3 

r ^ ^ 

> 

4- + 

In 

p 

► 

1 " 

a J 

a J 

J 


..2.31 


where p is the internal pressure arising from the misfit straining 
of the precipitate and cr^ is. the yield stress of the matrix. The 
model is applicable for both positive and negative stress-free 
misfit strain (c) values. These results indicate that the degree 
of plastic relaxation is strongly dependent on the precipitate 
size. The plastic relaxation is seen to increase with decrease in 
the relative stiffness of the precipitate and decrease with 
increase in the stress-free misfit strain. The plastic energy- 
contribution of the elasto-plastic deformation is shown to 
increase with increase in the relative stiffness of the 
precipitate . 

In the case of plastically deforming system following 
the linear strain hardening law, the total strain energy per unit 
volume of precipitate is given by [24] 


I - Kr 0-l)e^ - .2.32 


where, m is the linear strain hardening parameter i a s 1 1 c ^ 
and 0 = m/(l-m) , and f = 1 + 2 (l-v) . Other symbols have the 
same meanings as above . 
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In the case of systems exhibiting power law strain 
hardening, an iterative scheme has been presented to determine the 
stress and strain relations and hence the strain energy. The 
strain hardening law Act = k E)e” has been used, where n is 

plastic 

the strain hardening exponent (O^sn^sl) and k is a material 
constant . 

2.1.3 Finite Element models 

Brian et al [15] applied the finite element method to 

evaluate the elastic and plastic strain energies developed in a 
system due the presence of an ellipsoidal precipitate. The 
pioneering work in the application of finite element method to the 
evaluation of stress and strain energies during solid state phase 
transformation is done by Sen et al [13,14] . They developed a 

model to calculate the elastic [13] and elasto-plastic [14] strain 
energies of the matrix and the precipitate as a function of 

progress of transformation, where the shape of the precipitate was 
spherical. Precipitation from both the centre outwards and from 
the surface inwards were considered by them. 

2.2 STRAIN ENERGY AND HYSTERESIS IN METAL-HYDROGEN SYSTEMS 

Phase transformations in metal -hydrogen systems are 

manifested by charging hydride forming metals like niobium (Nb) , 
vanadium (V) with hydrogen, with accompanying formation of a 
hydride once the solubility limit of hydrogen in the metal is 
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exceeded. Reversal of the conditions result in dehydriding of the 
system. This process of hydriding and dehydriding shows a loss of 
efficiency in the form of hysteresis [15] . 

There three common modes of hysteresis in the 
metal -hydrogen systems [25] . They are: 

(i) Pressure hysteresis : Here the pressure needed for 
hydride formation, P^, is greater than that of hydride 
decomposition, . 

(ii) Solvus hysteresis : The terminal solid 
solubility of hydrogen in the host alloy is greater for hydride 
formation than for hydride decomposition. 

(iii) Thermal hysteresis : At a given solvus 
composition, the hydride-formation temperature is lower than the 
hydride-decomposition temperature . 

Extensive review of the subject has been done by 
Balasubramaniam [1] and is presented here. 

The origin of hysteresis in metal - hydrogen systems has 
been the topic of many studies starting as early as 1925 when 
Lambert and Gates [2 6] observed the phenomenon in the Pd-H system. 
Since then, several theories have been proposed to explain 
hysteresis in metal -hydrogen systems. Theses theories have been 
critically reviewed by Qian and Northwood [25] and the salient 
features of these theories would be explained here, briefly. Most 
of the theories emphasize the contribution of transformation 
strain as the possible origin of hysteresis. Ubbelhode [22] 
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proposed an additional strain factor in the phase rule to account 
for hysteresis. Scholtus and Hall [28] argued that the volume 
change on hydride formation would lead to a constraining stress on 
the matrix and they evaluated this stress to be approximately 
equal the yield stress of the matrix from which the hydride 
precipitates. Kuijpers and van Mai [29] proposed an explanation 
similar to that of Scholtus and Hall. The notion that plastic 
accommodation during hydride formation could be the possible 
source of hysteresis was proposed by Birnbaum et al [30] and 
extended by Flanagan and co-workers [31] . They proposed that the 
energy lost during the hysteresis cycle was accommodated by 
irreversible plastic deformation of the parent phase. Flanagan et 
al. have assumed that the plastic accommodation energy in the 
parent phase during hydride formation and that in the hydride 
phase during the hydride decomposition were the same. This is not 
a valid assumption because the elastic properties (and hence the 
dislocation formation energies) are significantly different for 
the metallic and hydride phases. Moreover, Birnbaum et al. also 
assumed that the elastic, chemical and surface accommodation 
energies were the same for hydride formation and decomposition 
and hence considered these energies to be reversible. While this 
may be true for chemical energy, it is not necessarily true 
correct for surface and elastic energies. Therefore, several 
simplifying arguments have been used by these workers while 
attributing the hysteresis effect to irreversible plastic 
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deformation. The other theories that were proposed to explain 
hysteresis are based upon saturation properties, chemical 
properties of the phases, crystallite size, metastable hydride 
formation, atomistic strain effect, elastic strain effects, and 
energy dissipation effects during hydrogen transformation. 
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CHAPTER 3 

FINITE ELEMENT MODEL 


3 . 1 INTRODUCTION 

The finite element method is a versatile tool for 
analyzing complicated engineering problems. In recent years it has 
been applied to the field of phase transformation in materials 
[11-15]. Using finite element method stress, strain and elastic 
and plastic accommodation energies developed during phase 
transformation can be predicted more accurately. 

This chapter briefly outlines the theory of Finite 
Element displacement formulation as applied to solution of 
problem of statics and the modeling of solid state phase 
transformation using finite element method. The two dimensional 
Finite Element displacement formulation of an elastic continuum is 
explained. Method of solution of plasticity problem using 
elasto-plastic analysis is stated in short . At the end Finite 
Element formulation of solid state phase transformation is 
explained . 

3.2 FINITE ELEMENT METHOD 

The analysis of an engineering system requires the 
idealization of the system into a form that can be analyzed, the 
formulation of governing equilibrium equations of this idealized 
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system, the solution of the equilibrium equations and the 
interpretation of the result. In many situations an adequate 
engineering model of the system is obtained by replacing the 
system by a finite number of well defined components. In contrast 
to this are the continuous systems, in which formulation leads to 
differential equations or equivalent statements, which imply an 
infinite number of elements. The exact solution of differential 
equations satisfying all boundary conditions is only possible for 
relatively simple systems, and numerical procedures must in 
general be employed to predict system response. In essence these 
procedures reduces the continuous system to a discrete 
idealization that can be analyzed in the same manner as discrete 
physical system. FEM is the most versatile numerical procedure 
for the approximation of continuum' problem. 

The Finite Element Method [32] essentially involves, 
dividing the domain into discrete or finite number of elements . 
The governing equations are reduced into simpler but approximate 
forms by discretisation. The resulting equations are replaced by 
suitable matrices to carry out the solution procedure. 

3.2.1 discretisation procedure 

The solution of stress and strain distribution in 
elastic continuum is the major field of application of Finite 
Element Method. In the Finite Element approximation the continuum 
is separated by imaginary lines or surfaces into a number of 
'finite elements'. The elements are assumed to be interconnected 
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at discrete number of nodal points situated on their boundaries . 
The displacement of these nodal points will be the basic unknown 
parameters of the problem. A set of functions is chosen to define 
uniquely the state of displacement within each 'finite element' in 
terms of its nodal displacements .The displacement functions define 
uniquely the state of strain within an element in terms of nodal 
displacements. These strains together with any initial strain and 
the constitutive properties of the material will define the state 
of stress throughout the element. 

For an elastic continuum, the governing differential 
equations can be obtained by minimizing the total potential energy 
of the system. The total potential energy , rr, can be expressed 
as [33] 


1 

n = 

2 

[cr] ^ [ e ] dv - 

[U]^[P] dv - 


V ^ 

V 


[U] [q] ds 


. .3.1 


Where 


O' -stress vector . 
e -strain vector. 

U -displacement at any point, 
p -body forces per unit volume, 
q -applied surface tractions. 

Integrations are taken over the volume V of the system 
and loaded surface area S. 
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The displacement U at any point inside an element is 
described as a function of nodal values 

U = [N] [U®] ' ..3.2 

Where [N] is a set of interpolation function termed as shape 
functions and [U^] is the vector of nodal displacement. 

The strain within the element is expressed in terms of 
nodal displacement as 

[e] =[B] [U®] ..3.3 

where [B] is the strain-displacement matrix, generally 
composed of derivatives of the shape functions. Stress is related 
to strains by use of an elastic matrix (property matrix) [D] , as 

[O'] = [D] [e] ..3.4 

When the shape function is such that no singularity exist in 
integrands of functional, the total potential energy of continuum 
will be the sum of the energy contributions of the individual 
elements . 


TT = S TT ..3.5 

e e 

where n is the total potential energy of element e. 

e 

Using eq (3 . 1) 
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TT 


U^]^[B]^[D] [B] [U®] dv - 


[U®]^ [N] ^ [p] 


dv - 


[U®] [N] [q] dS ..3.6 

■'s 

where Ve is the element volume and Se is the load surface area. 


Performance of the minimization of potential energy for element 
e with respect to nodal displacement results in 


Stt 

= [K®] [U®] - [I^ ] . .3.7 

au® 

where [F®] is the equivalent nodal forces for the element, 

[F®] = S [N]Up] dv + X [N]Uq] dS ..3.8 

V ^ V 

and [K] is the element stiffness matrix, 

i 

[K®] = X [B]Ud] [B] dv . .3.9 

, V 

The summation of the above terms in eq(3.7) over all the elements, 
when equated to zero results in a system of equilibrium equations 
for the complete continuum. These equations are then solved by any 
of standard technique to yield the nodal displacement . 

For an n noded element the shape functions is 
{N} ={ N^,N^ Nj 


..3.10 



The strain displacement relationships for an element in plane 
stress formulation are given by 


. o ^ ^ „ gu . dv 

X ax ' y ay ' 'xy“ ay ax 


..3.11 


which when expressed in matrix form gives equation 3.3, where 
strain displacement matrix is 



"n 

1 , X 

0 

N 

2 , X 


1 — 1 

td 
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1 > y 
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1 , X 
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2 , y 
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n , y 
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n , y n , X 


..3.12 


where u^ , v^ are the displacements in the global coordinates of 
the i^^ node of the element and 


N 

1 


X 


aNi 

ax 


N 


aNi 

1 , y sy 


..3.13 


In the axisymmetric formulation, the strain matrix is 



..3.14 


where e =(5u/5r),e = {5v/6z) , e _= (u/r) and r = [ (5u/5z) + (5V/5r) ] , 

r 2 o r z ^ 

with u and v the displacements in the polar coordinates and r the 
positional vector of the volume element of cross sectional area 
dA. 

The strain displacement matrix in axisymmetric formulation is 
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[B] 


N 

qI . r- 

N /r 
N 
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0 

N 

0 

. . . . N 
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N 

q 2 , r 
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. . . . o"’ 

N 

qI . 2 
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N 
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N 

1 , r 

3 , z 

2, r 
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n , r 


..3.15 


The property matrix for the plane stress formulation is 
E 


[D] = 




1 

V 

0 


V 

1 

0 


0 

0 


,3.16 


(l-i^) /2 

In the case of axisymmetric formulation, the property 


matrix is 


[D] = 


E (1-v) 


(l+i^) (l-2y) 


1 (y/(l-v)) (y/(l-y)) 0 

(y/(l-y)) 1 iv/il-v)) 0 

(y/(l-y)) (y/(l-y)) l 0 

0 00 (l-2y) / (2 (l-y) ) 


..3.17 

where, E is the Young's Modulus and v is the Poisson's ratio 
of the element material . 


3,2.2 plasticity formulation 


The difference in elastic and plastic behavior under 
uniaxial stress is closely brought out in Figure 4. In linear and 
nonlinear elastic behavior stress can be expressed as a function 
of strain. The main difference of plasticity formulation from 
nonlinear elastic formulation is that a unique stress- strain 
relationship is not available. Plastic deformation is essentially 
irreversible on unloading and incompressible in nature. 
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c .strain hardening plasticity 


d . Ideal plasticity 


Figure 4. Idealized elastic and plastic behaviour of a material under uniaxial stress. 
The up and down arrows indicate loading and unloading behaviours, respectively. 




The theory 

for 

plasticity 

formulation 

has been 

dealt 

with extensively by 

Owen 

and Hinton 

[33] . 

The 

theory has 

been 

presented below. 







In order 

to 

formulate 

a 

theory 

which 

models 


elasto-plastic material deformation, three requirements have to be 
met : 


1. An explicit relationship between stress and strain 
must be formulated to describe material behavior under elastic 
conditions, i.e. before the onset of plastic yielding. 

2 . A yield criterion indicating the stress level at 
which plastic flow commences must be postulated. 

3 . A relationship between stress and strain must be 
developed for post-yield behavior, i.e. when the deformation is 
made up of both elastic and plastic components . 

Before the onset of plastic yielding the relationship 
between stress and strain is given by the standard linear elastic 
expression: 


cr 

i J 


X 


i J k 1 


e 


k 1 


. .3.18 


where and are the stress and strain components 
respectively and is the tensor of elastic constants which 
for an isotropic material has the explicit form 


i j k 1 


= X 8 


1 j 


k 1 


+ (1 5 


1 k 


J 1 
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J k 


..3.19 


34 



where A and ju are Lame constants and 6 

1 J 

defined by 


is the Kronecker delta 


6 

i j 



i =j 

if 




..3.20 


The Von Mises yield criterion which has been used in the 

modeling, can be stated as, yielding occurs when , the second 

deviatoric stress invariant reaches a critical value, 

= k{k) ..3.21 

where k is a material parameter. 

/ 

J can be defined as 

2 

lrj2 .2 ,2, 2 2 2 

J = — [cr' + cr" + a' ]+ X + r + x ..3.22 

2 2 X y z xy yz xz 


where cr^, cr^ and are the principal deviatoric stresses in ' the 
X, y and z directions respectively in the cartesian coordinates, 
while X , X and x are the shear components . 

X y y z X z ■*" 


The yield criteria may be expressed as 

= •/3 k ..3.23 

where a- is the effective stress, given by 


O' 


/3 (j;) 


..3.24 


If a material deforms plastically, in general, the 
deformation would be composed of an elastic component and a 
plastic component. If the deformation load is applied 
incrementally and (dc^^)^ represents the plastic component of the 
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strain increment in the n^^ increment, then the increment in 
effective plastic strain is given by 

de = vTiTTT { (dc ) (de ) ..3.25 

p 1 j p i J p 

where the effective plastic strain increment is related to the 

effective stress by the relationship for uniaxial tension 

a- = a + dc ..3.26 

n n “ 1 p 

where H' denotes the strain hardening parameter (slope of the 
plastic portion of the stress-strain curve at the level of the 
effective stress at the n-l'^^ load increment) . 

After initial yielding, the material behavior will be 
partly elastic and partly plastic. During any increment of stress, 
the changes of strain can be divided into elastic and plastic 
components, so that 

de, = (dc, ) + (de, ,) ..3.27 

i J i J e 1 J P 

The elastic strain increment is related to the stress increment by 


dm' 


(de ) 

1 J e 


1 J 


(l-2l^) 


6 dm 

1 J kk 


..3.29 


24 E 

where E and v are the elastic modulus and Poisson' s ratio of 
the material, respectively, while dm' ^ and dm^^ are the deviatoric 
and the hydrostatic stress components respectively. 

The plastic component of the strain is formulated as 


follows. Taking f as the yield function, 

5f 


(de, ,) = dX 

i j P 


..3.30 


5m 


1 J 
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termed as the flow rule or normality condition since is a 
vector directed normal to the yield surface at the stress point 
under consideration. dA is the plastic multiplier (a 
proportionality constant) . Combining equations 3.28, 3.29 and 
3.30, we get 


(d€ 
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..3.31 


3.3 MODEL FOR TRANSFORMATION THROUGH PREDESTINED PATH 

In this section, modeling of the solid state transformation 
using finite element method is explained. The plane stress case 
has been taken to explain the principles. In the plane stress 
case, the non- zero stress and strain components are the two 
principal x and y components and one xy shear component on the x-y 
plane . 

3.3.1 Linear elastic analysis 

Stresses develop in the system during transformations 
due to the forces at matrix-precipitate interface. As the 
precipitate tries to occupy a space that is smaller (or larger) 
than the matrix from which it forms, the surrounding matrix (i.e. 
the untransformed region) resists this. Consider a metal -hydrogen 
system of finite size which is stress free at the external 
surfaces . 
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Assume that a region of matrix (radius r) is out out and 
allowed to transform freely to the hydride (Figure 5) . This free 
expansion leads to straining of the hydride. The developed strain 
is assumed to be constant in all directions and in the finite 
element analysis treated as a misfit or initial strain (e ) . In 

O 

the present case of the Nb^NbH transformation, the value of the 
misfit strain can be determined from the relation 


ppt 


matrix 


matrix 


..3.32 


where V is the molar volume of the precipitate and V is 

ppt matrix 

the molar volume of the matrix. The size of the precipitate would, 
in that case, increases in proportion to e . The resulting size of 


the precipitate would 

hence be 

o 

equal to 

[r(l+e^)] . 

If 

this 

transformed phase (of a 

larger size) is to be 

put back 

into 

the 

initial empty space 

of radius 

r in the 

matrix. 

a 

force 


proportional to would have to be applied to the precipitate. 
When the transformed phase is inserted back into the space from 
where it had been originally removed, the precipitates would apply 
a force on the surrounding matrix creating a strain field. The 
matrix and the precipitate tend to restrict their own elastic 
deformations, since both the matrix and the precipitate have 
non- zero elastic modulii. Therefore, at equilibrium, the 
precipitate would attain a size between r and [r(l+c )]. The 
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Figure 5. Definition of strain parameters used in the initial strain approach during 
forward hydrogen transformation for centre to surface path. 
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equilibrium position of the boundaries of each phase is indicated 
by dashed lines in Figure 5. It is interesting to note that the 
initial strain approach resembles the classic Eshelby approach 
[2,3] in determining the elastic strain energies of a misfitting 
precipitate. The forces due to the misfit strain act internally at 
the precipitate-matrix interface. Their nature are opposite during 
the forward and reverse transformations, due to the reverse nature 
of the misfit strains. Figure 6 shows the nature of the 
misfit- strain forces at the precipitate-matrix interface during 
forward and reverse hydrogen transformations . The system is 
divided into small elements for the purpose of finite element 
analysis. The equivalent mechanical force of each element, {F }, 

O 

required to create is given by [32] 



B D I dV 


..3.33 


where [B] is the strain-displacement matrix, [D] the material 
property matrix, {e } the initial -strain vector and V the volume 

O 

of the element. The superscript, T, refers to the transpose. The 
integration is performed over each element and the resulting force 
matrices are assembled to obtain the global load vector for the 
entire system. The strain-displacement matrix relates the nodal 
displacement vector {U} of each element to its total strain vector 
{e>. The total strain vector for each element, {c}, is obtained 
from the relation 
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..3.34 


{ . } = [b] { u } 


The nodal displacements and the initial -strain loads are related 
by 



. .3.35 


where [K] 
elemental 
formed by 


is the assembled stiffness matrix of the 

stiffness matrix t^lglement given in 

the equation 


system. The 
the matrix 




K 


element 


B ^ D B dV 


..3.36 


The unknown nodal displacement vector {U} is obtained by 
solving equation (3.35). The actual strain of each element is 
equal to the difference of the derived strain vector {e}, obtained 
from equation (3.34), and the initial -strain vector {c }. The 

O 

initial-strain vectors are zero for the untransformed elements, 
while they are non- zero for the elements representing the 
transformed region. Therefore, in general, the elemental stress 
vectors are obtained from the expression 



element 



element 




element 


..3.37 


where {a} is the transformation stress vector. 
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The elastic strain energy of the system ^^elastic^ 
then given by the relationship 
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elastic 
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r f 1 f 1 1 


cr 


e ^ e v 
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dV 


..3.38 


The total energy of the system is obtained by summing up the 
energies of individual elements . 


^no. of elements 

^elastic 2” ? _ 

1=1 

where the matrices have the same meanings as explained above but 
are now referenced to the elements. 

3.3.2 Elasto-plastic analysis 

The FEM formulation of elasto-plastic analysis using the 
incremental deformation theory is discussed below. The main 
advantage. of adopting the incremental deformation theory is that 
no iterative procedure is required in this case. The disadvantage 
is that very small load increments should be used. The matrix 
formulations used for the elasto-plastic analysis are similar to 
those used for the elastic axisymmetric formulation in the 
previous section, except that elasto-plastic stiffness have been 
used instead of the purely elastic stiffness used earlier. The 
initial-strain load is applied incrementally within each step of 
transformation . 

Incremental application of loading is essential for calculation 




V. 

1 


3 .39 
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of stresses and strains during plastic deformation. In the 
incremental loading procedure, the load matrix for the n^^ step of 
transformation is calculated only for the elements in the n*"^ 
layer of the mesh. This load is applied in p increments. The 
initial strain used in calculating the strain energies in every 
increment of loading (in each step of transformation) is 1/p*^*' of 
the initial strain calculated for the entire transformation step. 
The stresses and strains developed incrementally during a 
particular step of transformation are added to the previous ones, 
and the procedure continued. 

In incremental loading, the problem of material change 
need special attention, when the initial-strain approach has to be 
applied to a new layer of elements which earlier had the 
properties of the matrix but are now assigned the properties of 
the precipitate. The effective stresses, effective strains and 
total plastic strains for the the elements to undergo 
transformation are initialized before applying the loads because a 
new layer of elements undergoing transformation is assumed to be 
elastic at the beginning of the transformation step. In case the 
plastic deformation has occurred in the previous loading step, the 
accumulated plastic strain (up to the previous step of 
transformation) is added to the specified misfit strains to 
calculate the initial- strain load for the current step of 
transformation. The problem of material change during 
transformation is illustrated in Figure 7 which illustrates the 
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stress- strain curves used for Nb and NbH in the present analysis. 
NbH has been assumed to have a considerably high strain hardening 
coefficient such that it does not undergo yielding at any stage of 
phase transformation. Taking the strains to be instantaneously the 
same at the time the properties of the transforming elements 
change from that of the matrix to that of the precipitate, the 
stresses vary locally. During the forward transformation, a stress 
level at in the matrix would momentarily change to a in the 

precipitate, at the existing strain level, assuming elastic 
properties of the two phases (Figure 7) . Similarly, a stress level 
at a in the matrix would instantaneously change to b in the 
precipitate during the reverse transformation. The difference 
between the two stresses are converted into load vectors for the 
transforming elements using the following relation [14] 


■^Ao'l=AD-|e +e l 

l J L J I elastic plastlcj 


and 




dv 
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In situation where the initial-strain-load has to be 
applied in an element that was earlier defined as belonging to the 
matrix and had undergone plastic deformation, the plastic strain 
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total uniaxial strain 


e 



Figure 7. Schematic stress-strain curves of Nb and NbH showing the effect ot 
materal change. 
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must be subtracted, since this is the amount of strain already 
undergone by the element and it can only expand through balance of 
the total initial strain when the transformation occurs . The 
initial strain matrix is modified to account for the plastic 
deformation undergone by the element to transform. For 
axisymmetric formulation 


J 
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e 

- e 

f ] 
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1 } H 

0° 

- • 



rz{p) 


^ 0 

2 (p) > 


and for plane stress formulation 
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where e is the 

0 

initial strain calculated 

as (1/3) of 

the 

volumetric strain 

for 

free transformation 

of the matrix 

to 

precipitate and 

^r<p)' 

^0(p)^ ^rz(p) ' ^z(p) ' 

^xc p ) ' ^y (p) 

and 


) 9-^® the components of the total plastic strain undergone 
by the elements to undergo transformation in the n-1*^^ step. The 
initial load vector {F^} given by equation (3.33) for purely 
elastic behavior of the element is modified using equations 3.41 
and 3.43 to account for the plastic deformation and change of the 
element material from the matrix to the precipitate to give 
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The stresses for the transforming elements currently 
stored are updated to 





D 

precipitate 





The complete FEM solution procedure for the elasto-plastic case is 
given below. (The matrices have the same meanings as defined 
earlier. All symbols with A preceding it represent increment of 
property) : 

Step 1 : Calculation of total initial strain matrix {e } and 

O 

initial load matrix {Fq} for the transforming step using 
equation (3.45) . 

Step 2 : Application of part ( [AFq] ) of the load [F^] in each 
increment. Calculation of elemental stiffness with [D] replaced by 
[D , ] for the elements using 

elasto-plastic 

[ ]ela.a„t=| [ [ '> ] [ H 

Step 3 : Evaluation of displacements with the current load 
increment using 

[ 

where {AU} is the incremental displacement matrix for the current 
loading step. 




,3.48 
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step 4 : Calculation of incremental 


(elemental) stresses 


assuming elastic behavior by using 

{ } = [ “elastic ] [ { At } 
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where 


I Ac J- = B j I AU 


..3.50 


and {Act} is the incremental stress matrix. 

Step 5 : In the case of elements that have neither yielded 
earlier nor in this increment, incremental stresses calculated in 
step 4 are retained. For elements that have yielded in earlier 
increments but unloaded in the present one, the change in stresses 
in the present increment are taken as calculated in step 4 . For 
elements that have yielded or further plastically deformed during 
the current load increment, the stresses calculated in step 4 are 
reduced to the yield surface iteratively [33] . 

Step 6 ; The incremental stresses are added to previously 
stored ones (from earlier steps) . 

Step 7 : The cumulative elastic strains are obtained by the 
matrix relation 





'I 


^elastic ) , * 

element 
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“elastic 
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Step 8 : The incremental and total effective plastic strains 
c and stresses cr are calculated from the stress and strain 

p 

matrices obtained by step 7 from Von Mises criterion using [35] 
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^ = [i/(3/2)] {o-;^ 

and 

e = W{2/2)] {{c') {£' ) ..3.53 

P i J P i j P '' 


where (Cjj)p is the total plastic strain component and cr-'^ ^ is the 
deviatoric stress component. 

Step 9 : The cumulative elastic energy is obtained using 


'elastic 
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element 


element '■ 

while the incremental plastic strain energy is calculated from the 
following relation and added to the previously stored cumulative 
plastic strain energy 


no. of elements 

AE , . . = Z 

plastic ^ ^ 


O' Ae ..3.55 
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where cr is the effective stress of the element and 

effect Ive 

Ae is the change in effective plastic strain of the element 

plastic 

due to the load increment . 
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Yielding occurs if the effective stresses at a point exceeds the 
yield strength of the material at that point, where yield strength 
is given by 
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O' = cr° + e 
y y p 

where a° is the uniaxial yield strength of the material and E' is 
the strain hardening parameter. 

Step 10 ; Steps 2 to 9 are repeated till all the increments in each 
step have been applied. This is repeated for all the steps till 
completion of the transformation. 
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CHAPTER 4 

RESULTS AND DISCUSSION 


This chapter presents the accommodation energies 
obtained by the finite element model of hydrogen transformations 
presented in chapter 3. The results are presented in two sections. 
In the first section, the results obtained by linear elastic 
analysis are presented. The accommodation energies of the system 
obtained by plane stress formulation where the precipitate is in 
the shape of a thin disk and the accommodation energies obtained 
by axisymmetric formulation where the precipitate is a prolate 
spheroid are presented in that section. In the second section, the 
results of elasto-plastic analysis are presented. In this case 
also, energies of the system were calculated for disk shaped and 
prolate spheroid precipitates . Figure 8 shows the shape of the 
precipitate along with the element used to model the system in 
plane stress and axisymmetric analyses. The term aspect ratio 
refers to the ratio of semimajor axis to the semiminor axis in the 
case of disk shaped precipitate and in the context of 

axisymmetric analysis, the aspect ratio refers to the ratio a^ to 
a„ where a,, a..,, a., are the principle axes of prolate spheroid (a„ 
= a^) . A comparison of the accommodation energies for different 
aspect ratios (1, 2, 3 and 4) of the precipitate have been 
presented and discussed. 

The present model evaluates the accommodation energies 
during cc~^ hydrogen transformations m the Hb*“H and Pd“Ii (Appendix 
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Figure 8. Shape of precipitate and the element used to model the system in axisymmetric 
and plane stress analyses. The axes used are also indicated. 




A) systems by considering transformation progressing by two paths, 
viz precipitation starts at the centre and grows outwards towards 
the free surface (centre to surface path) and the precipitation 
starts at the surface of the system and grows inwards towards the 
centre (surface to centre path) . The a-/3 transformation in 
niobium-hydrogen system has been specifically addressed because 
the properties of the Nb and NbH and the experimental 
accommodation energies at 313 K are known [1] . (The results 
obtained for the hydrogen transformation in the palladium 
hydrogen system are given in Appendix A) . Niobium has a bcc 
structure (lattice parameter 3.3 A). The physical and mechanical 
properties of Nb and NbH used in the calculations are listed in 
Table 1 along with the reference source of data. The misfit volume 
strain was obtained from the molar volumes of Nb and NbH, and the 
linear misfit strain was assumed to be one third of the volumetric 
strain. The misfit strain during the forward transformation is 
positive and during 'the reverse transformation is negative. The 
misfit strain is found to be equal to 0.0439394 during a -> 0 

transformation and equal to -0.0420899 during 0 -» a transformation 
using equation [3.32]. The FEM calculations were performed on a 
HP-9000/850 computer. 

4.1 SELECTION OF FEM MODEL PARAMETERS 

While correct modelling of the physical problem into 
mathematical equations forms the core of the entire analysis, the 
s 0 X 0 c;tion of appropriate element and elemental mesh decides the 
accuracy of results. The judgment about the selection of various 
0 n\ 0 t;Q]^ j^ 00 cis 0 good deal of numerical experiment before a 



Table, 1. Physical and mechanical properties of Nb and NbH at 31 
K. Reference source is also indicated. 


Property 

Nb 

NbH 

Young's Modulus (E) 
(N/m^) 

1.22*10^^ 

[1] 

1.269*10^ ^ 

[1] 

Poisson's ratio (v) 

0.38 

0.38 


[1] 

[1] 

Yield stress (o-y) 

2.4828*10® 

4.3137*10® 

(N/m^) 

[1] 

[1] 

Crystal structure 

BCC 

Distorted 


[1] • 

BCC [1] 

Lattice parameter 

1 

3.3002 

3.445 

(A) 

[1] 

[1] 

Shear Modulus (fi) 

0.3752*10^^ 

0.452*10^ ^ 

(N/m^) 

[1] 

[1] 

Work hardening 

CO 

o 

H 

CD 

m 

100" • 

3.48*10 

parameter (H' ) 

(N/m^) 

[35] 


Molar volume { "v” ) 

10.8224*10"® 

12.3126*10"® 

(m^mol" ^ ) 

[1] 

[1] 


# Value assumed to prevent yielding in NbH 
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particular set of values is used for final analysis. In the 
present analysis, eight-noded isoparametric elements are used to 
model the system. A study of the sensitivity of the results with 
the type of element have been carried out using four-noded and 
eight-noded isoparametric elements and it was observed that both 
the type of elements yielded the same results . A similar 
observation was also noted by Sen et al [3 6] . The results obtained 
using different number of eight-noded isoparametric elements to 
model the system in the elastic axisymmetric analysis of forward 
transformation by the centre to surface path for aspect ratio 1 
are presented in Figure 9 . The calculated accommodation energies 
marginally depends upon the number of elements in to which the 
system is divided or, alternatively stated, the number of steps 
chosen to accomplish transformation. Calculations showed that 
using 10 steps of transformation or more, with the quadrant 
divided in to 4 radial sectors, gave converging results using 
either the eight-noded or the four-noded isoparametric elements. 

It has was earlier noted in Section 3.3.2 that in the 
elasto-plastic analysis, the load corresponding to the misfit 
transformation strain is applied in each stage of transformation 
in steps . The main advantage of adopting the incremental 
deformation theory is computer time- intensive iterative procedure 
is not necessary. In the incremental deformation theory used in 
the elasto-plastic analysis, the load should be applied in small 
increments so that in each load increment the amount of plastic 
deformation occurred is relatively small. The effect of varying 
the number of increments by which the load is applied during each 
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transformation stage on the system accommodation energy was 
studied. Figure 10 shows the system elasto-plastic accommodation 
energy during forward transformation by the center to surface path 
for different number of load increments (5, 10 and 20) used for 
the application of load in each stage of transformation. It is to 
be noted that the stress that is developed in each load increment 
for the 5 increment case is about 15% of the yield stress, which 
again is applied in a number of steps in each increment as 
explained by Owen et al [3.3] , while it is 10% and 3% of the yield 
stress when the number of increments are 10 and 20, respectively. 
There is no significant change in the energy of the system by 
changing the number of load increment steps. Therefore, in all the 
elasto-plastic analyses to be presented below, the number of load 
increments used in each transformation stage was 5 . 

In the present analyses the domain is divided into 
eighty elements and the transformation is modeled to occur in 2 0 
steps . Two types of domain used to model the system are shown in 
Figure 11. Figure 11a shows the cross sectional view of one 
quadrant of the ellipsoidal domain and Figure 11b shows the cross 
sectional view of one quadrant of the cylindrical domain. In the 
case of cylindrical system, the shape of the precipitate changes 
from that of the ellipsoid to cylinder on logarithmic scale with 
most of the shape change occurring in last few steps . Figure 12 
shows the results of elasto-plastic analysis of the forward 
transformation through the centre to surface path for the two 
different shapes of the system. The results obtained are similar. 
Therefore, the finite element mesh shown in Figure 11a is used to 
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Figure 10. Comparison of the system elasto-plastic accommodation energies during forward transformation 
through centre to surface path obtained by applying the load corresponding to the volume mismatch in different 
no of increments in each stage of transformation. 
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Figure 11. The finite element meshes along with the boundary conditions used for the calculation of accommodation 
energy in (a) an ellipsoidal domain and (b) cylindrical domain. 


ratio 1 
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Figure 12. Comparison of system elasto-plastic accommodation energies during forward transformation through 
the centre to surface path for different aspect ratios of an ellipsoidal precipitate in a cylindrical and an 
elhpsoidal matrix. 



model the system in the analysis to be presented henceforth. 


4.2 VALIDATION 

The results obtained by the linear elastic analysis of 
a-/3 transformation in Nb-H system using finite the element model 
have been compared with the accommodation energies obtained by 
applying other mathematical models of elastic accommodation. The 
models proposed by Eshelby [2] and Lee et al [23] are used for the 
comparison. In the model of Eshelby and Lee et al . , the matrix is 
an infinite continuum. Table 2 compares the elastic accommodation 
energies during during a-|3 hydrogen transformation in Nb-H system, 
wherein the precipitate is a spheroid. The reported FEM values are 
the elastic strain energies when the precipitate size tend to zero 
for the centre to surface path, the condition closely representing 
precipitation in an infinite continuum matrix. The magnitude of 
elastic accommodation energy in the Nb-H system depends upon the 
elastic strain energy model used. It was found that the results 
obtained by finite element method [14] is comparable with the 
results based on the models of Eshelby and Lee et al . 

4.3 LINEAR ELASTIC ANALYSIS 
4.3.1 Axisymmetric formulation 

The axisymmetric modeling of the ellipsoidal system was 
performed using eight -noded isoparametric elements (Figure 8) . 
Only one quadrant of the system has been considered because of 
symmetry. The domain is divided in to 20 concentric layers each 
consisting of four elements (Figure 11a) . The 20 concentric layer 
were transformed successively to depict 20 stages of the 



Table 2. Comparison of elastic accommodation energies during oc-p 
phase transformation in Nb - H system at 313 K. 


Model used for 


Elastic accommodation 

calculation. 

Formula 

energy (joules per mole 

[ ref ] 


metal in precipitate) 



forward 

reverse 

Finite Element 

! J-[<rl t {€) - dv 

4216 

4349 

Method 




Eshelby [2] 

24 (^) 

3a' V 

max 

3966 

4068 

Lee et al . 

[23] 

SUOi-yC Vmax 

a (y-1) +1 

3537 

4327 


^Definition of symbols used: 

4 = shear modules of matrix. 

V =molar valume . 

VV = |VNbH - VNb I 


y 


K 

pp t 


matrix 


K = Bulk Modules. 

V = poisson'B ratio 


a = 1+44 ^ / 3K ^ , 

matrix matrix 


a = 


(i+v ^ , ) 

matrix 

3(1-1^ ^ , ) 

matrix 


63 














transformation. All the elements in the transformed region were 
assigned the elastic properties (elastic modulii and poisson' s 
ratio) of the precipitate. The untransformed elements were 
assigned the elastic properties of the matrix. Calculation of the 
stresses, strains and strain energies were made at each stages of 
the transformation.. The load matrix at each stage of the 
transformation was modified according to equation [3.33] . At every 
stage, transformation was assumed to begin from the unstrained 
condition. Therefore, the energy calculated at every stage of 
transformation is the cumulative energy until that stage of the 
transformation. The transformation progressing from centre of the 
matrix outwards (centre to surface path) and from surface of the 
matrix inwards (surface to centre path) are addressed in the 
present analysis. The accommodation energies are calculated for 
different aspect ratios of the precipitate. 

As the matrix is stress free at the beginning and at the 
end of the transformation, the strain energies will be zero at 
beginning and end of the transformation irrespective of the 
aspect ratio. This is true for all transformation paths* At any 
intermediate stage of transformation, there will be two phases 
each constraining the expansion/contraction of other phase. 
Therefore, there will be a net strain energy at any intermediate 
stage of transformation. 

Figure 13 shows the variation of elastic accommodation 
energy of the system as a function of progress of transformation 
by centre to surface transformation path for transformation. 
The energy is expressed in joules per mole of Nb in the system. 



The advantage of using this energy unit has been discussed 
elsewhere [12] . Figures 14a and I4b show the elastic strain 
energies stored in the matrix and the precipitate, respectively, 
for various aspect ratios for the forward transformation. Since 
NbH is assumed to be a stochiometric monohydride, one atom of 
hydrogen atom is associated with an atom of the metal in the 
hydride. Therefore, the number of moles of Nb in the system is a 
constant for given domain. The abscissa is the mole fraction of 
NbH, which increases with progress of transformation in the 
transformation and decreases with progress of transformation in 
/3-4a transformation. Figure 14b shows that the precipitate energy 
increases with aspect ratio whereas the matrix energy (Figure 14a) 


and total 

accommodation energy 

(Figure 13) 

decrease 

with 

an 

increase 

in aspect ratio 

of 

the 

precipitate . 

For 

the 

transformation progressing 

from 

the 

centre 

to surface. 

the 


precipitate NbH will be constrained by a surrounding shell of 
matrix. When the aspect ratio is one, i.e. precipitate is a 
sphere, the state of stress inside the precipitate is purely 
hydrostatic. In this case the precipitate constraint on the matrix 
is high, there by making the elastic deformation of matrix 
relatively large compared to the elastic deformation of 
precipitate. Hence, the matrix energy is dominant in this case. 
For aspect ratios other than one, the precipitate is no longer in 
a state of pure hydrostatic stress. There will be shear stresses 
present inside the precipitate which increase with increase in 
aspect ratio. Therefore, the constraint offered by the precipitate 
decrease with increasing aspect ratio, which results in a 
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Figure 13. Variation of system accommodation energy as a function of fraction transformed during forward 
transformation through centre to surface path for different aspect ratios of the precipitate using elastic 
axisymmetric analysis. 
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Figure 14. Variation of accommodation energy stored in the (a) matrix and 
(b) precipitate during forward transformation through centre to surface 






decrease in the energy of matrix with increasing aspect ratio 
(Figure 14a) . The strain energy in the precipitate increases 
with an increase in aspect ratio since elastic deformation inside 
the precipitate will be more in this case (Figure 14b) . In order 
to understand the effect of shear components on the accommodation 
energies, the average effective stress in the precipitate and the 
matrix at 50% transformations for the forward and reverse 
transformations are tabulated in Table 3 . The effective stress cr 
is defined as 


a = V3 (J ) 
2 


1/2 


4.1 


where J is the 
2 

j' = - 
2 2 


second deviatoric stress invariant. 


[cr' ' 

X 


+ cr 


,2 


+ cr' 
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+ T 


xy 


+ T 


y z 


+ X 


.4.2 


Effective stress provides information about the existence of shear 
components of stress. It can be seen from Table 3, that the 
effective stresses in the precipitate is negligible, when the 
aspect ratio of the precipitate is one and increases as the aspect 
ratio increase, whereas the effective stresses inside the matrix 
decreases as the aspect ratio increases. For the centre to 
surface transformation path, the major contribution to the total 
accommodation energy is from the matrix, so the total 
accommodation energy of the system decreases with increase in 
aspect ratio. As the precipitate approaches the external boundary- 
stresses get relaxed. Therefore, the strain energy decreases 
during the last stages of transformation for all the aspect ratios 
(figures 13 through 16) . Eventually, the strain energy decreases 



Table 3. Average effective stress (GPa) inside the precipitate and 
the matrix at 50s transformation for different aspect ratios of the 

precipitate . 




_ Average Effetctive stress (GPa) 



precipitate 

matrix 



aspect ratio 

12 3 4 

aspect ratio 

1 2 3 4 - 

centre to 
surface 

forward 

reverse 

~ 0 1.50 1.68 1.71 

~ 0 1.46 1.62 1.70 

6.31 5.54 5.34 5.29 

6.29 5.50 5.29 5.20 

surface to 
centre 

forward 

5.98 5.14 4.93 4.88 

- 0 1.75 1.98 2.04 


reverse 

5.51 4.77 4.57 4.52 

~ 0 1.67 1.90 1.96 
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to zero as all of the specimen is transfoimaed and there are no 
more misfit constraints present in the system. 

For the reverse transformation starting from centre of 
the system, the variation of the accommodation energies stored in 
the system, the matrix and the precipitate, which are depicted in 
figures 15, 16a and 16b, respectively, show the same pattern as 
in the forward transformation. The above results shows that for 
centre to surface transformation path, whether it is forward or 
reverse transformation, most of the elastic strain energy is 
accommodated in the outer phase surrounding the inner phase. 

The parabolic relationship of the strain energy for the 
case of spheroid precipitate (figure 13 through 20) is not 
surprising and it is expected regardless of the relaxation. If w 
is the elastic strain energy (per unit volume of precipitate) and 
(j) is the volume fraction of the precipitate phase, then the total 
strain energy is given by oxp {!-<()) , which is the well known 
parabolic relationship [37] . Strictly speaking cj is a function of 
<p, but it . is relatively insignificant as the transformation 
strain/stress field is of short range. 

Figure 17 shows the variation of system accommodation 
energy with progress of transformation during forward 
transformation through the surface to centre path, while Figures 
18a and 18b provide the matrix and the precipitate energies, 
respectively. In the case of transformation starting from surface 
growing inwards, a shell of precipitate material grows into the 
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Figure 15. Variation of system accommodation energy during reverse transformation through the centre to surface 
path for different aspect ratios of the precipitate using elastic axisymmetric analysis. 
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(b) precipitate energy 

Figure 16. Variation of accommodation energy stored in the (a) matrix and 
(b) precipitate during reverse transformation through centre to surface 
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Figure 17. Variation of system accommodation energy during forward transformation through the surface to centre 
path for different aspect ratios of the precipitate using elastic axisymmetric analysis. 
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(a) matrix energy 
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Figure 18. Variation of accommodation energy stored in the (a) matrix and 
(b) precipitate during forward transformation through surface to centre path 




matrix. The results show that the energy of the matrix (Figure 
18a) increases with increase in the aspect ratio of the 

precipitate, whereas the precipitate energy (Figure 18b) and the 
total energy of the system (Figure 17) decreases as the aspect 
ratio of the precipitate increase. When the aspect ratio is one, 
the matrix will be under pure hydrostatic state of stress. In 
this case the contribution from the matrix to the total energy of 
the system is small. As the aspect ratio increases, shear 

stresses develop inside the matrix thereby increasing its elastic 
deformation (Table 3) and, therefore. The matrix energy increases 
with aspect ratio (Figure 18a) . The constraint on the precipitate 
by matrix decreases with an increase in the aspect ratio since it 
is no longer in the hydrostatic state of stress. This reduces the 
strain energy inside the precipitate shell (Figure 18b) , which is 
the major contributor to the total energy. The total energy 
decreases with an increase in the aspect ratio of the precipitate 
(Figure 17) . 

, The accommodation energy of the system during reverse 

transformation through surface to centre path is shown in Figure 
19 and the accommodation energies stored in the matrix and the 
precipitate are shown in Figure 20a and Figure 20b, respectively. 
The energies show the same variation as in the forward 
transformation through surface to centre path. As in the case of 
centre to surface transformation, most of the strain energy in the 
surface to centre case is accommodated in outer phase surrounding 
the inner phase (i.e. NbH in the forward transformation and Nb in 
the reverse transformation) . 
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Figure 19. Variation of system accommodation energy during reverse transformation through the surface to centre 
path for different aspect ratios of the precipitate using elastic axisymmetric analysis. 




ACCOMMODATION ENERGY (J/mol Nb in the system) ACCOMMODATION ENERGY (J/mol Nb in the systen 



MOL FRACTION NbH 


(a) matrix energy 



MOL FRACTION NbH 

(b) precipitate energy 

Figure 20. Variation of accommodation energy stored in the (a) matrix and 
(b) precipitate during reverse transformation through surface to centre path 




4.3.2. Plane stress formulation 

Eight -noded isoparametric elements have been utilized to 

model the system in the case of plane stress fprmulation. The 

system is specified with a thickness of IjLim and a cross sectional 

-3 2 

area of 7.8539x10 m (which corresponds to the area of a 
circular system of radius 0.05m). The precipitate is disk shaped. 
Only one quadrant of the system has been considered due to 
symmetry. Each quadrant has been divided in to 80 elements with 
281 nodal points. 

The calculated variation of the elastic strain energy as 
a function of fraction of matrix that has been transformed to 
precipitate for the forward (oc->i3) and reverse O-^a) hydrogen 
transformations by the centre to surface and surface to center 
paths are presented in Figures 21 and 22, respectively. The 
energy is expressed as joules per mole of metal in the system. The 
abscissa denotes the volume fraction of hydride in the system at 
each stage of transformation. 

The elastic strain energies shows the same variation as 
that in axisymmetric case. However, as the constraints to the free 
expansion is lesser in plane stress case compared to the 
axisymmetric case, the numerical values of the strain energy at 
any particular stage of transformation is lower in plane stress 
case when compared to axisymmetric case. 

4.4 ELASTO-PLASTIC ANALYSIS 

The finite element model described in the earlier 
section to evaluate the elastic strain energies during room 
temperature a-/3 hydrogen transformation yields results comparable 
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(a) forward transformation 



(b) reverse transformation 

Figure 21. Variation of system accommodation energy during (a) forward and 
(b) reverse transformation through centre to surface path for different aspect 
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(b) reverse transformation 


Figure 22. Variation of system accommodation energy during (a) forward and 
(b) reverse transformation through surface to centre path for different aspect 





to the elastic strain energies obtained with other mathematical 
models [2,17] . However, the calculated strain energy in the Nb-H 
system would be largely reduced due to plastic deformation which 
is expected to occur due to the relatively large misfit 
strain (~4%) involved during the a-|3 transformations in the Nb-H 
system [1] . Transmission microscopic studies of hydride formation 
in niobium [31,38] have revealed that a phase deforms plastically 
during both the forward and reverse hydrogen transformation. The 
effect of plastic deformation in lowering elastic strain energy 
during solid phase transformation has been earlier studied by Lee 
et al by considering the matrix to be ideally plastic [23] and to 
strain harden according to the power law [24] . However, these 
models cannot be used to estimate the accommodation energy if the 
precipitate deforms plastically as is the situation during the 
reverse hydrogen transformation through centre to surface 

path in metal hydrogen systems. 

The present FEM model has been utilized to evaluate the 
elasto-plastic accommodation energies during the 
transformations in the Nb-H system. The results obtained by both 
the plane stress analysis (where the precipitate is disk shaped) 
and the axisymmetric analysis (where the precipitate is an ob 
spheroid) are presented below. The effect of changing the aspect 
ratio on the accommodation energy of the system is specific y 
addressed. NbH is assumed to be elastic as it has been 
experimentally observed that |3 phase does not yield ' g 
hydrogen transformation in Nb-H system [39] . 



4.4.1. Axisymmetric formulation 

The axisymmetric modelling of the ellipsoidal (prolate 
spheroid) system was performed using eight-noded isoparametric 
elements. The dimensions of the prolate spheroid were chosen such 
that the total volume is same for all aspect ratios. The system is 
specified by a volume of 5.235e-4 m^ which corresponds to a sphere 
of radius 0.05m. It has been observed that a small change in 
volume does not affect the accommodation energies which have been 
normalized by dividing the energy by the volume of the system. The 
evaluation of the accommodation energies have been performed for 
different aspect ratios (l, 2, 3 and 4). Only one quadrant of the 
system has been considered due to the symmetry. Each quadrant was 
divided in to 80 elements with 281 nodal points. The 
transformation was modelled to occur in 20 steps. In each step, 
the load corresponding to the volume mismatch due to 
transformation was applied in 5 steps. 

In the elasto-plastic analyses, NbH is assumed to be 
elastic in all the transformations, and hence the plastic 
deformation,, if any, is confined to only Nb, The results of the 
elasto-plastic finite element analysis- are presented in Figures 
23 through 32. The strain energies are expressed as a function 
fraction of matrix transformed. Figure 23 compares the 
accommodation energies for the forward transformation through 
centre to surface path obtained by elastic and elasto-plastic 
analysis in the case of a spherical precipitate. The total 
accommodation energy in the elasto-plastic case is lower than the 
slastic strain energies generated in the elastic analysis except 
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Figure 23. Comparison of system accommodation energy obtained by elastic and elasto-plastic axisyrametric 
analyses during forward transformation through centre to the centre to surface path. 




at the end of transformation. This reduction in the total energy 
in the elasto-plastic case is due to the plastic deformation in 
the matrix. The energy m the matrix due to plastic deformation 
is irreversible which is expended during the transformation. This 
appears as the residual energy at the end of the transformation in 
the elasto-plastic case. The present results are in conformity 
with that of Sen et al [14] . 

The system energy obtained by the elasto-plastic 
analysis of the forward transformation assuming centre to surface 
path is 9^^® Figure 24, while the matrix and precipitate 

energies are provided in Figures 25a and 25b. Here NbH, which is 
assumed to be elastic, precipitates inside the matrix (Nb) and 
grows outwards . The total strain energy of the system for 
different aspect ratios of the ellipsoidal precipitate is shown in 
Figure 24 which is the sum of plastic energy of the matrix and 
elastic energy of both matrix and precipitate. The energy of the 
matrix and precipitate are shown in Figures 25a and 25b, 
respectively. The major contribution to the system strain energy, 
at any stage of transformation, is from the matrix. Figure 25b 
shows an increase in accommodation energy of the precipitate with 
an increase in aspect ratio. When the aspect ratio is one the 
precipitate is under hydrostatic state of stress. As the aspect 
ratio is increased, shear stresses also develop inside the 
precipitate and therefore, strain energy of the matrix is 
composed of elastic and plastic component, of which plastic strain 
energy is dominant . The computed values of deviatoric components 
of stress inside matrix, which causes the plastic deformation 
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Figure.24 Variation of system elasto-plastic accommodation energy during forward transformation through centre 
to surface path for different aspect ratios of the precipitate using axisymmetric formulation. 




(a) matrix energy 



(b) precipitate energy 

Figure 25. Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during forward transformation through centre to 
surface path for different aspect ratios of the precipitate using axisymmetric 



the strain 


increases as the aspect ratio increases. Therefore, 
energy of the matrix increases with increase in aspect ratio 
(Figure. 25a) . The strain energy of the matrix is large relative to 
that of the precipitate, and hence the total energy of the system 
increases as aspect ratio increases (Figure 24) 

The variation of the strain energy in the case of 
precipitation of Nb (reverse transformation) from the centre of 
the system outwards is presented in Figures 26 through 28. Here 
the matrix phase (NbH) is assumed to be elastic. The energies are 
given as a function of progress of transformation and it is to be 
noted that at the beginning of transformation the mole fraction of 
the NbH (abscissa) is one and it decreases as the transformation 
progresses and finally reaches zero when the whole matrix is 
transformed to Nb. The matrix strain energy (Figure 27a) which is 
purely elastic, decreases with increase in aspect ratio in the 
first few stages of transformation. This variation is similar to 
the results obtained by the linear elastic analysis (Figure 16a) . 
The strain energy of the precipitate (Nb) is the sum of elastic 
and plastic strain energies of which the elastic part is 
significant. The elastic and plastic part of the strain energy of 
the precipitate is shown in Figure 28. The state of stress is not 
uniform about the 0-axis and this nonuniformity becomes more 
pronounced as the aspect ratio increases. It has been observed 
that the plastic part of the strain energy stored in the 
precipitate increases with increase in aspect ratio, whereas the 
elastic strain energy decreases. This decrease in elastic part of 
the total energy of the phase in the initial stages of 
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Figiire.26 Variation of system elasto-plastic accommodation energy during reverse transformation through the 
centre to surface path for different aspect ratios of the precipitate using axisynunetric foimulation. 




ACCOMMODATION ENERGY (J/mol Nb in the s; 



(a) matrix energy 



(b) precipitate energy 

Figure. 27 Variation of elastO"plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during reverse transformation through the centre 
to surface path for different aspect ratios of precipitate using axisymmetric 




(a) elastic strain energy 



(b) plastic strain energy 

Figure 28. Variation of (a) elastic and (b) plastic part of accommodation 
energy stored in the Nb precipitate during reverse transformation through 
the centre to surface path for different aspect ratios of the precipitate using 



transformation with increasing aspect ratio is found in all the 
situations where the phase is plastically deforming. The elastic 
energy is maximum for the precipitate having aspect ratio one and 
decreases as the aspect ratio increases . The above picture is true 
only for the initial stages of transformation. The non-uniform 
plastic deformation of the precipitate leads to the locking up of 
elastic strains inside the precipitate. The effect of locked up 
elastic strains is more pronounced in the latter stages of 
transformation, i.e when considerable plastic deformation has 
already occurred. Since the stress distribution will be more 
uniform as aspect ratio is nearly one, the locked up elastic 
strains will be smaller for these cases. In the latter stages of 
transformation, the elastic strain energy of the precipitate is 
small in the case of aspect ratio equal to one and increases as 
aspect ratio increase (Figure 28a) . The plastic part of strain 
energy of the precipitate increases with an increase in aspect 
ratio as shown in Figure 2 8b. The variation of the total energy of 
the precipitate (Figure 27b) is same as that of elastic part of 


the strain 

energy . 

The locked 

up elastic 

strain 

inside 

the 

precipitate 

which 

is present 

in the 

latter 

stages 

of 


transformation offer constraint to the matrix. The effect of this 
constraint is to increase the elastic strain energy of the 
matrix for higher aspect ratios in the latter stages of 
transformation as shown in Figure 27a. Figure 26 shows that the 
pattern of variation of total accommodation energy of the system 
is same as that of the NbH matrix. But at the end of 
transformation there will be net energy lost due to the 



irreversible plastic deformation. 

The elasto-plastic analysis of the transformation 
has also been performed when the transformation begins from the 
surface and the precipitate grows inwards (surface to centre 
path) . The results of this analysis for the forward transformation 
is given in Figures 29 and 30. In this case NbH precipitates as 
a hollow sphere. The precipitate will not yield at any stage of 
transformation due to its high work hardening parameter and high 
yield strength. When the precipitate is a hollow sphere (aspect 
ratio 1), the matrix is under hydrostatic state of stress. 
Therefore there will be no plastic deformation of the matrix in 
this case. As the aspect ratio increases, the strain energy inside 
the precipitate, which is purely elastic, decreases (Figure 30b) 
as obtained in the linear elastic analysis (Figure 18b) . The 
strain energy of the matrix (Nb) for aspect ratio other than one 
has two components, elastic strain energy and plastic strain 
energy. It has been observed that the plastic part of strain 
energy of the matrix increases with an increase in aspect ratio 
whereas the . elastic strain energy of the matrix decreases with an 
increase in aspect ratio. Note that this is pattern of variation 
is same as the variation of precipitate energy that obtained for 
initial stages of the reverse hydrogen transformation from centre 
to surface. The variation of the total energy of the matrix is 
shown in Figure 30a. The total accommodation energy of the system 
which is the sum of the energies of matrix and the precipitate is 
shown in Figure 29. The total energy follows the path of 
precipitate energy in the initial stages of transformation but at 
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Fit^ure.29 Variation of system elasto-plastic accommodation energy during forward transformation through the 
surface to centre path for different aspect ratios of the precipitate using axisymmetric formulation. 



ACCOMMODATION ENERGY (J/mol Nb in the system) ACCOMMODATION ENERGY {J/mol Nb in the S] 



(a) matrix energy 



(b) precipitate energy 

Figure 30. Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during forward transformation through surface to 
centre path for different aspect ratios of the precipitate using axisymmetric 

fnriTnilation. 



the end of transformation total energy is equal to plastic part of 
the strain energy . The results shows that for aspect ratio other 
than one there is not much change in the strain energy at the end 
of transformation with change in aspect ratio. 

Figures 31 and 32 present the accommodation energies 
during the reverse transformation (NbH->Nb) progressing from the 
surface towards centre of the system. The precipitate in this case 
is Nb which yields due to the volume misfit strain as in the case 
of forward transformation through centre to surface path. Figure 
32b shows the variation of the precipitate (Nb) energy during the 
transformation. The elastic energy of the precipitate is found to 
be negligibly small in this case. As the aspect ratio increases 
the plastic part of the strain energy inside the Nb precipitate 
shell increases and therefore the total strain energy of the 
precipitate increases. The inner NbH matrix here is assumed to be 
elastic throughout the transformation. The energy of the matrix 
(NbH) is shown in the Figure 32a which shows that the energy of 
the matrix increases with the aspect ratio similar to the results 
obtained by linear elastic analysis of the reverse transformation 
through surface to centre path (Figure 20a) . The total energy of 
the system is presented in Figure 31 as a function of fraction 
transformed and the variation is similar to the variation of 
energies of the precipitate shell which is prime contributor to 
the total energy. The strain energy at the end of transformation 
is found to be more or less same for all the aspect ratios other 
than one . 

Table 4 presents the residual energies present at the 
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Figure.31 Variation of system elasto-plastic accommodation energy during reverse transformation through the 
surface to centre path for different aspect ratios of the precipitate using axisymmetric formulation. 




(a) matrix energy 



(b) precipitate energy 

Figure.32 Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during reverse transformation through surface to 
centre path for different aspect ratios of the precipitate using axisymmetric 

fnrmnintion 



Table 4. Elasto-plastic accommodation energy (joules per mol 
metal ) at the end of hydrogen transformation in Nb-H system 
obtained by axisymmetric formulation. 


Transformation 

path 

aspect 

ratio of 

the precipitate 

1 

2 

3 

4 

centre to surface: 

forward 

214.71 

224.64 

232.22 

235.40 

reverse 

0.00 

315.97 

446.73 

469.07 

surface to centre: 

forward 

0.00 

73.00 

70.04 

78 .47 

reverse 

437.25 

439.06 

437.98 

446.56 
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end of transformation for all the cases described above. In the 
case of a spherical precipitate (for aspect ratio 1) , the strain 
energy at the end of transformation are zero for reverse 
transfoinmation through centre to surface path and forward 
transformation through surface to centre path. In both cases, 

niobium is under hydrostatic state of stress. The maximum energies 
are obtained for the reverse transformations wherein the final 
phase is plastically deformed niobium. It can also be noticed in 
the table that the accommodation energies increase with increasing 
aspect ratio except in the case of reverse transformation by the 

surface to centre to centre path. This could be due to the 

presence of shear components in increasing magnitude as the aspect 
ratio is increased which results in a larger amount of plastic 
deformation energies with increasing aspect ratios. The 

accommodation energy at any given aspect ratio in the case forward 
transformation is significantly lower for the surface to centre 
path when compared to the centre to surface path. As noticed 
earlier, in the case of elastic analysis, most of the energy is 
accommodated by the outer phase that surrounds the inner phase due 
to the greater deformation suffered by it. When the deformable 
metallic matrix (Nb) surround the inner phase (as in the case of 
forward transformation through the centre to surface path) , it 
plastically deforms to a greater extent than when the metallic 
matrix is present as the inner phase (as in the case of forward 
transformation by the surface to centre path) . Therefore, the 
energies are down for the latter case when compared to the former. 
The Table 4 provides the accommodation energies at the 
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end of transformation. In order to understand the behaviour at 
midway in the progress of transformation, the individual 
accommodation energies of the precipitate and matrix phases for 
the centre to surface forward and surface to centre reverse 
transformations at 50% transformation are presented in the Table 
5. Note that, at 50% transformation, in both these cases, the 
inner phase NbH is enveloped by the outer phase Nb. The results in 
Table 5 indicates that, as already noted, the system 
accommodation energy in the elastic case decreases with increasing 
aspect ratio and this due to the major contribution of the outer 
phase (Nb) , whose accommodation energy shows a similar trend in 
the variation with aspect ratio. As already noted, this is 
reflected in the decrease in average effective stress in the 
outer (Nb) phase with increasing aspect ratio (Table 3) . In the 
elasto-plastic analysis, it is first observed that the 
accommodation energies are lowered when compared to the elastic 
case and this is due to plastic deformation, as has been shown 
analytically by Lee et al [23] . However, here the total 
accommodation energy increases with increase in the aspect ratio 
which is essentially due to a similar variation in the energy of 
the outer phase (Nb) with aspect ratio. Most of the accommodation 
energy of the outer phase (Nb) is plastic in nature. The increased 
plastic accommodation energy in the Nb matrix when it occurs as 
the outer phase is due to increased yielding of this phase with 
increase in aspect ratio due to the generation of additional shear 
stresses in this phase. This leads to the total accommodation 
energy increasing with increasing aspect ratio as plasticity 



Table 5. Accommodation energies (Joules per mol Nb) stored in the 
matrix and the precipitate at 50% transformation for different 
aspect ratios of the precipitate, where the inner phase is NbH and 
outer phase is Nb. 
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161 

165 
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115 

138 
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916 

832 

789 

771 
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840 

763 

724 

707 
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11 
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189 

197 

203 

206 
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forward 
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124 

123 

126 

127 
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determines the accommodation energies in the elasto-plastic case. 


4.4.2 plane stress formulation 

The thickness of the specimen used for the experimental 
investigation of reversible hydrogen transformations, for example 
by the pressure-composition-temperature (PCT) technique, is very 
small compared to its other dimensions. The transformation in such 
a system can be modelled as a plane stress problem, which reduces 
the complexity of the model while maintaining sufficient accuracy. 
The results obtained by the plane stress analysis of hydrogen 
transformation in the Nb-H system are presented below. Here the 
the precipitate is disk shaped. Eight-noded isoparametric elements 

have been used to model the system. The system is specified by a 
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thickness of 1 jjxa and a cross sectional area of 7.8539x10 m. Due 

to symmetry only one quarter of the system has been considered, 

which is divided in to 20 layers each having 4 elements. The 

loading procedure is similar to that used in axisymmetric 

analysis, i.e. five incremental steps for each stage of 

transformation . 

The system accommodation energy during the forward 
transformation progressing from centre of the system outwards is 
presented in Figure 33a. Note that the total energy of the system 
increases with an increase in the aspect ratio. The variation of 
accommodation energy with progress of transformation and the 
variation of energy with respect to aspect ratio are similar to 
the results obtained in elasto-plastic analysis of the 
axisymmetric ellipsoidal system. The strain energy at the end of 



transformation in this case is smaller compared to the residual 
energy obtained by axisymmetric analysis. 


Figure 33b depicts the variation of the accommodation 
energy during reverse hydrogen transformation through the centre 
to surface path. In the axisymmetric case, the precipitate was 
under hydrostatic stress for aspect ratio one, whereas here the 
precipitate has shear component even for aspect ratio one. 
Therefore, plastic deformation occurs in the precipitate for all 
the aspect ratio. The results shows that among the four aspect 
ratio considered, the maximum residual energy is obtained for 
aspect ratio 3 . Apart from this the variation of the total energy 
with aspect ratio follows the same pattern as is obtained in 
axisymmetric analysis. 

In the case of forward transformation through surface to 
centre path, the dependence of the accommodation energy on aspect 
ratio shows the same trend (Figure 34a) as that was obtained for 
an ellipsoidal precipitate. However, the matrix energy in this 
case will have plastic component for aspect ratio one since it is 
not under pure hydrostatic state of stress as in the case in 
axisymmetric analysis. For aspect ratio other than one, there is 
no significant change in the residual energy of the system with a 
change in aspect ratio . It is to be noted that most of the energy 
is accommodated in the matrix (Mb) as seen in Figure 35. The 
plastic and elastic components of the accommodation in the Nb 
matrix is provided in Figure 36 which shows that most of the 

energy in the Nb matrix is plastic in nature. 

The total strain energy of the system during the reverse 
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(a) forward 



MOL FRACTION NbH 

(b) reverse 


Figure 33. Variation of system elasto-plastic accommodation energy during 
(a) forward and (b) reverse transformations through the centre to surface path 
for different aspect ratios of the precipitate using plane stress formulation. 




(a) forward 



(b) reverse 

Figure 34 Variation of system elasto-plastic accommodation energy during 
(a) forward and (b) reverse transformations through the surface centre path 
for different aspect ratios of the precipitate using plane stress formulation. 





(a) matrix energy 



(b) precipitate energy 

Figure 35. Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during forward transformation through surface to 
centre path for different aspect ratios of the precipitate using plane stress 
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(a) elastic strain energy 
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(b) plastic strain energy 

Figure 36. Variation of (a) elastic and (b) plastic part of acconomodation 
energy stored in the Nb matrix during forward transformation through the 
surface to centre path using plane stress formulation. 



transformation by the surface to centre path increases with 
increase in aspect ratio at the initial stages of transformation 
(Figure 34b) . However, at the end of transformation the strain 
energies are almost same for aspect ratio other than one which is 
in conformity with the results obtained in axisymmetric analysis. 
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CHAPTER 5 
SUMMARY 


5.1 CONCLUSIONS 


The salient conclusions of the present finite element 
study are the following. The initial strain approach of the finite 
element method (FEM) has been successfully applied in order to 
evaluate the effect of precipitate aspect ratio on the 
accommodation energies during solid state hydrogen transformations 
in the Nb-H system. The transformation was modelled to occur from 
the surface of the system inwards and also from the centre of the 
system outwards . One of the main advantages of the FEM model is 
that finite systems can be modelled. The FEM method was utilized 
to evaluate the elastic and elasto-plastic accommodation energies 
during hydrogen transformations for the plane stress and 
axisymmetric conditions. The individual (elastic and plastic) 
energies stored in the matrix and precipitate phases can be 
evaluated by finite element analysis (FEA) . The effect of the 
aspect ratio of the precipitate on the accommodation energy has 
been primarily addressed. 

In the elastic FEA, it was observed that the 
accommodation energy of the system decreases as the aspect ratio 
of the precipitate increases, using both the plane stress and 
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axisymmetric formulations . This was found to be true for both the 
forward and reverse transformations by the center to surface and 
surface to centre transformation paths. In all these cases, the 
accommodation energy stored in the inner phase increases with 
increasing aspect ratio while the accommodation energy stored in 
the outer phase decreases with increasing aspect ratio. 
Therefore, the elastic analysis reveals that most of the energy is 
accommodated in the phase which surrounds the inner phase. 

In the elasto-plastic FEA, the accommodation energies 
during the course of the transformation are generally lowered 
compared to that obtained by the elastic analysis. This is due to 
plastic deformation that occurs in the metallic Nb phase. Residual 
energies are present at the end of transformation which are 
essentially due to the irreversible plastic deformation that 
occurs during the transformation. Exception to the above are when 
the metallic phase is the inner phase having aspect ratio one in 
which case the metallic phase does not yield as the state of 
stress is purely hydrostatic. The accommodation energies at the 
end of transformation are higher in the case of the 
transformations in which the final phase is plastically deforming 
niobium compared to the transformations in which the final phase 
is niobium hydride. The accommodation energies generally increase 
with increasing aspect ratio for all the transformation paths 
considered. 

The finite element model has also been applied to obtain 
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the elasto-plastic accommodation energies for hydrogen 
transformations in the Pd-H system. The nature of the results are 
similar to that obtained for the Nb“H system. 

5.2 SCOPE FOR FUTURE STUDY 

The following aspects could be considered in detail while 
^PPlyiJ^9’ the finite element method to evaluate accommodation 
energies during solid state transformations: 

a) The model applied in the present thesis assumes a single 
precipitate growing in a single crystal. This could be extended to 
take into consideration multiple nucleation and growth events 
(i.e. interacting precipitates) and also a multigrained matrix. In 
this case, grain boundary effects would also have to be considered 
and suitably accounted for in the FEM model. 

b) The transformation is assumed to progress through a 
pre-destined transformation path whereas the growth of the 
precipitate, in practice, would be governed by the stress state in 
the system, i.e. the precipitate would prefer to grow in the 
direction where the accommodation energy density is favorable. 
This should be incorporated into the model. 

c) The transformations could be modelled in three dimensions by 
choosing suitable three dimensional elements (for example, brick 
or tetrahedral elements) . 

d) Several other shapes of the precipitate could be modelled by 


111 



the FEM technique . 

e) It was implicitly assumed that the misfit strain was one third 
of the volumetric strain for the transformations considered in the 
present work. However, the misfit strains would be different in 
different crystallographic directions due to specific orientation 
relationships between the hydride and matrix. This could be easily 
incorporated into the present model (when applying for the Zr-H 
system) as only the input parameters have to be suitably modified. 
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APPENDIX A 


FEM RESULTS FOR THE Pd-H SYSTEM 

Palladium/hydrogen system has been extensively 
experimentally investigated for the hysteresis behaviour [l] . In 
terms of practical use, hydrogen is most often associated with 
palladium because of the wide spread use of the metal as a 
membrane which permits a rapid penetration of hydrogen and 
provides a convenient method of separating it out of gaseous 
mixtures. From an academic stand point, interest in the system has 
centered primarily around the attempts to provide theoretical 
explanations in the form of pressure -composition relationship 
which shows the existence of a hysteresis between absorption and 
desorption p-c relationship [1] . This appendix presents the 
elasto-plastic accommodation energies obtained for the hydrogen 
transformations in palladium-hydrogen system evaluated using the 
Finite element method. 

The FEM model of hydrogen transformation described in 
chapter 3 is used to model the a-(3 transformation in the Pd-H 
system. The mechanical and physical properties of palladium and 
palladium hydride which were used for the calculation are given in 
Table Al. The work hardening parameter for the Pd was obtained 
from the stress-strain curve presented in Ref. [4] . The hardening 
is almost linear in plastic range. The work hardening parameter H' 
is obtained from the relation l / = 1/ E + 1 / H' . The misfit 
strain for the a->|3 transformation is 0.03506 and that for p->a 
transformation is -0.0317 [2] .elasto-plastic accommodation 
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energies during hydrogen transformation were evaluated assuming 
the precipitate to be a prolate spheroid. The domain was divided 
in to 80 isoparametric eight-noded elements. The transformation is 
modeled to take place in 20 steps. All the transformation paths 
considered in the case of Nb-H system which are described in 
chapter 4 were modelled here also. Palladium hydride is assumed to 
be elastic in all the cases. 

The results of the finite element analyses are shown in 
Figure A1 through Figure A6. The figures show the variation of the 
accommodation energies as function of the fraction transformed for 
different aspect ratios of the precipitate. Figure A1 shows the 
system accommodation energy for (a) forward and (b) reverse 
transformation through the centre to surface path, while Figure A2 
shows the accommodation energies stored in the matrix and the 
precipitate for the forward transformation and Figure A3 shows the 
accommodation energies stored in the matrix and precipitate for 
reverse transformation. The accommodation energies in the case of 
surface to centre path are shown in Figures A4 through A6. The 
variation of the accommodation energy with fraction transformed 
and with aspect ratio of the precipitate are same as that has been 
previously obtained in the case of Nb-H system, which were 
presented in chapter 4. However the numerical values of the 
energies are smaller due to lower value of misfit strain and also 
due to the difference in the physical and mechanical properties of 

the two systems . 
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ACCOMMODATION ENERGY (J/mol Nb m the system) ACCOMMODATION ENERGY (J/mol Nb in the s: 



MOL FRACTION NbH 

(a) forward transformation 



MOL FRACTION NbH 


(b) reverse transformation 


Figure Al. Variation of system elasto-plastic accommodation energy during 
(a) forward and (b) reverse transformations through centre to surface path for 
different aspect ratios of the precipitate using axisymmetric formulation. 





MOL FRACTION NbH 


(a) matrix energy 



MOL FRACTION NbH 


(b) precipitate energy 

Figure A2. Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during forward transformation through centre to 
surface path for different aspect ratios of the precipitate using axisymmetric 
ff^rmnlnlion. 



(a) matrix energy 



(b) precipitate energy 

Figure. A3 Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during reverse transformation through the centre 
to surface path for different aspect ratios of precipitate using axisymmetric 




(a) forward transformation 



(b) reverse transformation 

Figure A4. Variation of system elasto-plastic accommodation energy during 
(a) forward and (b) reverse transformations through surface to centre path for 
different aspect ratios of the precipitate using axisymmetric formulation. 



(a) matrix energy 



(b) precipitate energy 

Figure A5. Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during forwaid transformation through surface to 
centre path for different aspect ratios of the precipitate using axisymmetric 

-FrM-rm 1 1 Cl r\n 



in the system) 



Figure. A6 Variation of elasto-plastic accommodation energy stored in the 
(a) matrix and (b) precipitate during reverse transformation through surface to 
centre path for different aspect ratios of the precipitate using axisymmetric 






Table. A1 Physical and mechanical properties of Pd and PdH 
room temperature. Reference source is also indicated. 


Property 

Pd 

PdH 

Young's Modulus (E) 
(N/m^) 

1.25*10^^ 

[3] 

1.15*10^^ 

[3] 

Poisson's ratio (i^) 

0.39 

0.39 


[3] 

[3] 

Yield stress (vy) 

0.52*10® 

0.52*10^®° 

(N/m^) 

[5] 


Crystal structure 

FCC 

FCC 


[2] 

[2] 

Lattice parameter 

3.894 

4.026 

(A) 

[2] 

[2] 

Work hardening 

0.12*10® 

100 

0.12*10 

parameter (H' ) 

(N/m^) 

[4] 


Molar volume ( V ) 

8.89079*10"® 

9.82594*10"® 

{m^mol~ ^ ) 

[2] 

[2] 


# Value assumed to prevent yielding in PdH 
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